/*
   FieldVis - A Tool for Visualizing Astrophysical Magnetohydrodynamic Flow
   
   Copyright (C) 2006 - 2008     FieldVis Development Team
  
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
  
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
  
   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

   File        : Field.cpp
   Author      : Jeremy Moc  (mocj@msi.umn.edu)

   Description : The vector and scalar field interface to the data sets.
*/

#include "Field.h"

///////////////////////////////////////////////////////////////////////////////

float ScalarField::GetInterpolatedValue(const float x, const float y, const float z) const
{
   int xFloor = (int) x;
   int yFloor = (int) y;
   int zFloor = (int) z;

   int iMax = (xFloor >= xDim - 1 ? 0 : 1);
   int jMax = (yFloor >= yDim - 1 ? 0 : 1);
   int kMax = (zFloor >= zDim - 1 ? 0 : 1);
   
   float value = 0.0;
   float weight;

   for (int i = 0; i <= iMax; i++)
      for (int j = 0; j <= jMax; j++)
         for (int k = 0; k <= kMax; k++)
         {
            weight = fabs( (1 - (x - xFloor + i)) *
                           (1 - (y - yFloor + j)) *
                           (1 - (z - zFloor + k)) );
            
            value += weight * GetGridPointValue(xFloor + i, yFloor + j, zFloor + k);
         }

   return value;
}


///////////////////////////////////////////////////////////////////////////////
Vector3D VectorField::GetInterpolatedVector(const float x, const float y, const float z) const
{
  int xFloor = (int) x;
  int yFloor = (int) y;
  int zFloor = (int) z;

  int iMax = (xFloor >= xDim - 1 ? 0 : 1);
  int jMax = (yFloor >= yDim - 1 ? 0 : 1);
  int kMax = (zFloor >= zDim - 1 ? 0 : 1);
   
  float xValue, yValue, zValue;
  float weight;

  xValue = yValue = zValue = 0.0;

  for (int i = 0; i <= iMax; i++)
     for (int j = 0; j <= jMax; j++)
        for (int k = 0; k <= kMax; k++)
        {
           weight = fabs( (1 - (x - xFloor + i)) *
                          (1 - (y - yFloor + j)) *
                          (1 - (z - zFloor + k)) );
            
           xValue += weight * xField->Get(xFloor + i, yFloor + j, zFloor + k);
           yValue += weight * yField->Get(xFloor + i, yFloor + j, zFloor + k);
           zValue += weight * zField->Get(xFloor + i, yFloor + j, zFloor + k);
        }

   return (Vector3D(xValue, yValue, zValue));
}
